Setting

path_proj = here::here()
path_source = file.path(path_proj, "source")

source(file.path(path_source, "simulation", "simulations_functions_final.R"))
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
## 
##     smiths
source(file.path(path_source, "functions", "plot_function.R"))
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
source(file.path(path_source, "functions", "fit_function.R"))
source(file.path(path_source, "functions", "table_function.R"))

# place for draws
# mac
posterior_draws_path = file.path(Sys.getenv("HOME"), "Desktop", "draws", "testEach")

#Windows
#posterior_draws_path = file.path(Sys.getenv("USERPROFILE"), "Desktop", "draws", "testEach")

#data path
data_save_path = file.path(path_proj, "data", "fitted_model", "simulation", "4. b_ou")
#models
q_constant <- file.path(path_proj, "source", "models", 
                      "q-constant.stan")
b_constant <- file.path(path_proj, "source", "models", 
                     "b-constant.stan")
b_rw <- file.path(path_proj, "source", "models", 
                     "b-rw1.stan")
b_ou <-  file.path(path_proj, "source", "models", 
                      "b-ou.stan")

compiled_models <- list(
  q_constant = cmdstan_model(q_constant),
  b_constant = cmdstan_model(b_constant),
  b_rw = cmdstan_model(b_rw),
  b_ou = cmdstan_model(b_ou)
)

models_to_use <- c("q_constant", "b_constant", "b_rw", "b_ou")
###### setting #####

# data
alpha_increase_seq_1 <- seq(10, 750, by = 30)
alpha_decrease_seq_1 <- seq(750, 10, by = -30)
alpha_lamb =  c( rep(10,5), alpha_increase_seq_1 + rnorm(alpha_increase_seq_1,10,10), 
                 alpha_decrease_seq_1 + rnorm(alpha_decrease_seq_1,10,10),
                 rep(10,5))
beta_lamb = 0.5
T = 60
# reprot delay
D <- 15;

# Time period for checking
D_check <- 5

first_date <- as.Date("2024-01-01")

scoreRange <- c(first_date+days(9), first_date+days(19), first_date+days(29),
                first_date+days(39), first_date+days(49))

Fully reported case

Simulation

params_b_ou_FR <- list(
  data = list(
    alpha_lamb = alpha_lamb,
    beta_lamb  = beta_lamb,
    T       = T,
    date_start = as.Date("2024-01-01"),
    D = D
  ),
  q_model = list(
    method        = "b_ou",
    method_params = list(theta_logb = 0.3, mu_logb = log(0.7), init_logb = log(0.7), sigma_logb = 0.2,
                         theta_logitphi = 0.3, mu_logitphi = 1, init_logitphi = 1, sigma_logitphi = 0.2)
  )
)
b_ou_FR <- simulateData(params_b_ou_FR)
## [1] "here we go"
par(mfrow = c(2, 1))
plot(b_ou_FR$b, pch = 19, type = "b")
plot(b_ou_FR$phi, pch = 19, type = "b")

par(mfrow = c(1, 1))
matplot(t(b_ou_FR$q), type = "l", lty = 1, ylim = c(0, 1))

### Exploratory analysis

# exploritary analysis
page_num <- ceiling(nrow(b_ou_FR$case_reported_cumulated)/16)
exp_plot_ou <- fit_exp_plot(b_ou_FR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# print(exp_plot_ou)
# 
# exp_plot_ou$coefficients
exp_b_data_ou<- data.frame( date = as.Date(rownames(b_ou_FR$case_reported_cumulated)),
                          b = exp_plot_ou$coefficients$b)

exp_b_plot_ou <- ggplot(exp_b_data_ou, aes(x = date, y = b)) +
  geom_point(color = "black", size = 1.5) +       
  geom_smooth(method = "loess", se = TRUE,        
              color = "blue", fill = "grey", alpha = 0.5) +
  theme_minimal() +
  labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")

print(exp_b_plot_ou)
## `geom_smooth()` using formula = 'y ~ x'

Model fitting

# out_b_ou_FR <-
#   nowcasting_moving_window(b_ou_FR$case_reported_cumulated, scoreRange = scoreRange,
#                           case_true = b_ou_FR$case_true,
#                           start_date = as.Date("2024-01-01"),
#                           D = D,
#                           methods =models_to_use,
#                           compiled_models = compiled_models,
#                           iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
#                           num_chains = 3, suppress_output = T,
#                           posterior_draws_path = file.path(posterior_draws_path, "b_ou")
#                           )
# 
# 
# save(out_b_ou_FR, file = file.path(data_save_path, "FR_b_ou.RData"))
load(file.path(data_save_path, "FR_b_ou.RData"))

results_b_ou_FR <- nowcasts_table(out_b_ou_FR, D = D, report_unit = "day", 
                          methods = models_to_use)

plots_b_ou_FR <- nowcasts_plot(results_b_ou_FR, D = D, report_unit = "day", methods = models_to_use)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
for (i in 1:length(results_b_ou_FR)) {
  print(calculate_metrics(results_b_ou_FR[[i]],
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 85.22 70.55 62.37 69.49         129.22           0.5 q_constant
## 2 19.56  7.57  9.24  5.64          51.91           1.0 b_constant
## 3 23.04  8.33 10.69  5.84          58.30           1.0       b_rw
## 4 19.91  7.60  9.32  5.45          57.60           1.0       b_ou
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 55.60 11.34 38.79 10.77          99.56          0.80 q_constant
## 2 61.74  7.69 25.93  4.09          77.86          0.85 b_constant
## 3 50.65  6.28 21.14  3.30          93.11          0.90       b_rw
## 4 59.77  7.35 23.78  3.63          91.81          0.95       b_ou
##     RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 207.01 14.54 71.48 6.73         112.68          0.83 q_constant
## 2 190.78 13.20 59.94 4.99         107.87          0.87 b_constant
## 3  24.01  1.84 10.57 1.18         141.98          1.00       b_rw
## 4  54.73  3.82 16.97 1.55         147.04          1.00       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 47.16  4.30 25.56 2.81         114.28          0.95 q_constant
## 2 39.89  3.60 19.18 2.06         111.69          0.95 b_constant
## 3 25.96  2.42  9.50 1.08         132.64          1.00       b_rw
## 4 26.00  2.38  8.43 1.00         131.17          1.00       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 17.38  2.37 12.72 1.83         107.65             1 q_constant
## 2 13.68  2.16  9.82 1.48         106.45             1 b_constant
## 3 12.17  3.01  4.85 1.06         111.39             1       b_rw
## 4  9.97  2.48  4.14 0.94         110.63             1       b_ou
for (i in 1:length(results_b_ou_FR)) {
  print(calculate_metrics(data.table::last(results_b_ou_FR[[i]],D_check),
                          methods = models_to_use))
}
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 119.21 69.53 108.52 68.11         222.22           0.2 q_constant
## 2  27.65 10.03  17.66  7.59          84.82           1.0 b_constant
## 3  32.57 11.38  20.69  8.63          96.81           1.0       b_rw
## 4  28.14 10.18  17.88  7.50          95.60           1.0       b_ou
##     RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1  91.84 12.78 79.53 11.12         195.01           0.6 q_constant
## 2 123.18 15.02 94.31 11.86         149.82           0.4 b_constant
## 3 101.17 12.37 78.88  9.89         207.62           0.6       b_rw
## 4 119.45 14.53 89.65 11.15         202.43           0.8       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 504.83 34.74 364.40 25.57         262.63           0.2 q_constant
## 2 466.25 32.02 319.70 22.28         246.81           0.2 b_constant
## 3  57.57  4.04  49.10  3.53         444.63           1.0       b_rw
## 4 133.58  9.17  89.63  6.27         472.82           1.0       b_ou
##     RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 115.26 10.27 90.21 8.12         199.41           0.6 q_constant
## 2  98.66  8.76 67.85 6.10         193.81           0.8 b_constant
## 3  71.91  6.56 52.79 4.85         350.24           1.0       b_rw
## 4  72.29  6.44 48.32 4.36         337.03           1.0       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 22.96  4.86 20.81 4.26         122.62             1 q_constant
## 2 22.97  5.39 20.94 4.52         119.81             1 b_constant
## 3 37.28  9.35 27.93 6.63         166.41             1       b_rw
## 4 30.57  7.61 24.48 5.68         158.41             1       b_ou
print(plots_b_ou_FR)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

### Parameter Checking

# try the third one, "2024-01-30"
T_test = T * 3/6

# b_ou
varnames_b_ou <- out_b_ou_FR$b_ou[[3]]$summary()$variable


mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("sigma_log_b"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("sigma_logit_phi"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("theta_log_b"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("theta_logit_phi"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("mu_log_b"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("mu_logit_phi"), prob = 0.95, prob_outer = 0.95)

param_true = tibble(
    parameter = grep("^b\\[.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$b[1:T_test]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("b"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^phi\\[.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$phi[1:T_test]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("phi"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^lambda\\[.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$lambda[1:T_test]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("lambda"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^q\\[10,.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$q[10,]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws(grep("^q\\[10,.+\\]$", varnames_b_ou, value = TRUE)), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^N\\[.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("N"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

Not fully reported case

Simulation

# ou_NFR
params_b_ou_NFR <- list(
  data = list(
    alpha_lamb = alpha_lamb,
    beta_lamb  = beta_lamb,
    T       = T,
    date_start = as.Date("2024-01-01"),
    D = D
  ),
  q_model = list(
    method        = "b_ou",
    method_params = list(theta_logb = 0.2, mu_logb = log(0.1), init_logb = log(0.1), sigma_logb = 0.15,
                         theta_logitphi = 0.2, mu_logitphi = 1.5, init_logitphi = 1.5, sigma_logitphi = 0.15)
  )
)

b_ou_NFR <- simulateData(params_b_ou_NFR)
## [1] "here we go"
par(mfrow = c(2, 1))
plot(b_ou_NFR$b, pch = 19, type = "b")
plot(b_ou_NFR$phi, pch = 19, type = "b")

par(mfrow = c(1, 1))
matplot(t(b_ou_NFR$q), type = "l", lty = 1, ylim = c(0, 1))

Exploratory analysis

# exploritary analysis
page_num <- ceiling(nrow(b_ou_NFR$case_reported_cumulated)/16)
exp_plot_ou <- fit_exp_plot(b_ou_NFR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
# print(exp_plot_ou)
# 
# exp_plot_ou$coefficients
exp_b_data_ou<- data.frame( date = as.Date(rownames(b_ou_NFR$case_reported_cumulated)),
                          b = exp_plot_ou$coefficients$b)

exp_b_plot_ou <- ggplot(exp_b_data_ou, aes(x = date, y = b)) +
  geom_point(color = "black", size = 1.5) +       
  geom_smooth(method = "loess", se = TRUE,        
              color = "blue", fill = "grey", alpha = 0.5) +
  theme_minimal() +
  labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")

print(exp_b_plot_ou)
## `geom_smooth()` using formula = 'y ~ x'

Model fitting

# out_b_ou_NFR <-
#   nowcasting_moving_window(b_ou_NFR$case_reported_cumulated, scoreRange = scoreRange,
#                           case_true = b_ou_NFR$case_true,
#                           start_date = as.Date("2024-01-01"),
#                           D = D, 
#                           methods =models_to_use,
#                           compiled_models = compiled_models,
#                           iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
#                           num_chains = 3, suppress_output = T,
#                           posterior_draws_path = file.path(posterior_draws_path, "b_ou")
#                           )
# 
# 
# save(out_b_ou_NFR, file = file.path(data_save_path, "NFR_b_ou.RData"))
load(file.path(data_save_path, "NFR_b_ou.RData"))

results_b_ou_NFR <- nowcasts_table(out_b_ou_NFR, D = D, report_unit = "day",
                          methods = models_to_use)

plots_b_ou_NFR <- nowcasts_plot(results_b_ou_NFR, D = D, report_unit = "day", methods = models_to_use)

for (i in 1:length(results_b_ou_NFR)) {
  print(calculate_metrics(results_b_ou_NFR[[i]],
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 12.54 27.74 10.36 22.44          83.51           1.0 q_constant
## 2 63.09 50.82 41.74 48.92          44.51           0.3 b_constant
## 3 65.19 50.87 42.71 49.00          47.51           0.3       b_rw
## 4 65.18 51.47 43.00 49.72          46.60           0.3       b_ou
##     RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1  97.53 20.11 59.79 16.87         105.21          0.75 q_constant
## 2 107.32 21.24 68.16 18.51          98.11          0.70 b_constant
## 3 101.82 19.10 61.49 15.82         151.86          0.85       b_rw
## 4  97.64 19.29 59.39 16.16         137.12          0.85       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 230.91 25.77 164.53 23.15          97.91          0.33 q_constant
## 2 151.09 17.19  91.79 13.17         112.18          0.60 b_constant
## 3 125.71 15.59  68.59 10.97         230.08          0.87       b_rw
## 4 149.25 16.58  84.42 11.93         181.58          0.83       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 208.06 22.57 157.32 19.46         105.31          0.32 q_constant
## 2 122.22 14.75  80.34 11.51         124.76          0.65 b_constant
## 3 120.17 14.57  64.99  9.33         246.94          0.95       b_rw
## 4 110.66 13.29  58.07  8.31         246.22          0.95       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 188.98 22.04 152.18 19.94          99.57          0.24 q_constant
## 2  93.19 12.37  66.32 10.09         114.71          0.60 b_constant
## 3  43.68 10.30  31.00  6.52         224.60          1.00       b_rw
## 4  42.72  9.40  30.72  6.10         211.42          1.00       b_ou
for (i in 1:length(results_b_ou_NFR)) {
  print(calculate_metrics(data.table::last(results_b_ou_NFR[[i]],D_check),
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 16.91 23.72 15.79 17.89         147.01           1.0 q_constant
## 2 88.46 41.99 72.15 40.41          76.00           0.4 b_constant
## 3 91.46 42.97 74.21 41.21          81.82           0.4       b_rw
## 4 91.43 43.46 74.62 41.89          80.00           0.4       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 187.03 23.63 178.25 22.74         210.24           0.2 q_constant
## 2 203.78 25.76 195.49 24.98         195.03           0.2 b_constant
## 3 195.01 24.79 186.64 23.92         350.04           0.4       b_rw
## 4 187.88 23.68 179.15 22.82         303.84           0.4       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 346.27 24.98 334.78 24.28         214.63           0.0 q_constant
## 2 188.77 13.66 164.09 12.00         248.82           0.4 b_constant
## 3 159.53 11.67 129.36  9.58         735.84           1.0       b_rw
## 4 207.92 15.28 174.11 12.86         517.85           0.8       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 356.28 34.65 342.81 32.82         168.42           0.0 q_constant
## 2 237.49 23.66 202.05 19.77         201.00           0.2 b_constant
## 3 308.36 30.28 286.99 27.68         534.86           0.6       b_rw
## 4 288.14 28.41 262.57 25.42         469.82           0.6       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 125.73 25.42 124.44 25.22         116.81           0.0 q_constant
## 2  49.77 10.41  46.79  9.69         139.01           0.8 b_constant
## 3  79.11 16.99  75.93 15.97         386.05           1.0       b_rw
## 4  82.07 17.20  79.81 16.53         276.44           1.0       b_ou
print(plots_b_ou_NFR)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Parameter Checking

# try the third one, "2024-01-30"
T_test = T * 3/6

# b_ou
varnames_b_ou <- out_b_ou_NFR$b_ou[[3]]$summary()$variable


mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("sigma_log_b"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("sigma_logit_phi"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("theta_log_b"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("theta_logit_phi"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("mu_log_b"), prob = 0.95, prob_outer = 0.95)

mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("mu_logit_phi"), prob = 0.95, prob_outer = 0.95)

param_true = tibble(
    parameter = grep("^b\\[.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_NFR$b[1:T_test]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("b"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^phi\\[.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$phi[1:T_test]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("phi"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^lambda\\[.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$lambda[1:T_test]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("lambda"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^q\\[10,.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$q[10,]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws(grep("^q\\[10,.+\\]$", varnames_b_ou, value = TRUE)), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^N\\[.+\\]$", varnames_b_ou, value = TRUE),
    x = b_ou_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("N"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)